*Install packages
*ssc install randcmd
*ssc install texsave

* Set path

if "`c(username)'"=="x" { /
	global dir =  "C:/Users//EcoTurismo/"
	global base_out "C:/Users/EcoTurismo\CreatedData\"
	global outcomes "C:/Users/EcoTurismo/Paper/Tables/"
	global graphs "C:/Users/EcoTurismo/Paper/Images/"
	global numbers "C:/Users//EcoTurismo/Paper/Numbers/"
}

* Import data
use "${base_out}/cleaned_Ecoturismo_followup", clear


************ Table 1: Summary statistics *************************************** 
* Panel A: Municipalities
* output: Tables/sum_balanceMuni_1
preserve
duplicates drop codmpio, force
replace distanciaBog=0 if distanciaBog==-1 | distanciaBog==.
replace distanciaBog=distanciaBog*70

label var homicidios "Homicides 2018 per 100k inh"
label var distanciaBog "Distance to Bogotá (km)"
label var areaBosque2017k2 "Forest cover area (km2)"
label var nbi2005 "Poverty index"
label var mean_business "\% of business"

global tab_muni "areaBosque2017k2 nbi2005 homicidios mean_business"
	
texdoc init "$outcomes/sum_balanceMuni_1.tex" , replace force 
loc varN=1
cap mat drop coefs
foreach var of varlist $tab_muni {
	
	keep if time == 0
	
	local ytitle : var label `var'
				
		qui sum `var' if treat==1
		local mean1 : di  %4.2f r(mean)
		local sd1: di %4.2f r(sd)
		local N1: di %4.0f r(N)
		local N_total = r(N)
		
		qui sum `var' if treat==0
		local mean0 : di  %4.2f r(mean)
		local sd0: di %4.2f r(sd)  
		local N0: di %4.0f r(N)
		local N_total: di %4.0f `N_total' + r(N)
		
		ttest `var' , by(treat)
		local ttest: di %4.2f `r(p)'
		local diff: di %4.2f `mean1' - `mean0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
	
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\
	if mod(`varN',4)==0 {
	tex   \midrule
	tex         N                 &   `N1' &    `N0'       &    `N_total'     \\  
	tex   \midrule
	}
	loc varN=`varN'+1		
}	
texdoc close
restore

* Panel B: Business Baseline
* output: Tables/sum_balanceMuni_2                                              
preserve

keep if time==0
keep if in_sample_ntourist == 1

*label var ntourist "N tourists normal month"
label var nworkers "N workers"
label var LossSitio2Km "Defo area (ha) around site (2km)"
label var M1_201A "N local tourist"
label var M1_201B "N foreign tourist"

qui sum ntourist if treat==1
local N1: di %4.0f r(N)
local N_total = r(N)		
qui sum ntourist if treat==0
local N0: di %4.0f r(N)
local N_total: di %4.0f `N_total' + r(N)
		
global tab_muni "M1_201A M1_201B nworkers M1_205 LossSitio2Km" 
	
texdoc init "$outcomes/sum_balanceMuni_2.tex" , replace force 
loc varN=1
cap mat drop coefs
foreach var of varlist $tab_muni {
	
	local ytitle : var label `var'
				
		qui sum `var' if treat==1 
		local mean1 : di  %4.2f r(mean)
		local sd1: di %4.2f r(sd)
		*local N1: di %4.0f r(N)
		
		qui sum `var' if treat==0 
		local mean0 : di  %4.2f r(mean)
		local sd0: di %4.2f r(sd)  
		*local N0: di %4.0f r(N)
		
		
		ttest `var' , by(treat)
		local ttest: di %4.2f `r(p)'
		local diff: di %4.2f `mean1' - `mean0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
	
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\
	loc varN=`varN'+1		
}
	tex   \midrule
	tex         N                 &   `N1' &    `N0'       &   `N_total'   \\  
	tex   \midrule
texdoc close
restore

* Panel C: Hosehold Baseline
* output: Tables/sum_balanceMuni_3                                              
preserve

keep if in_sample_reg_ecot

keep if time==0

*label var M1_114D "Gender (Male)"
label var LossCasa2Km "Defo area (ha) around HQ (2km)"

	    qui sum M1_504A if treat==1 										
		local N1: di %4.0f r(N)		

		
		qui sum M1_504A if treat==0
		local N0: di %4.0f r(N)

global tab_muni "M1_504A month_income LossCasa2Km"						// La variable M1_114D es sexo y no todos pusieron
	
texdoc init "$outcomes/sum_balanceMuni_3.tex" , replace force 

loc varN=1
cap mat drop coefs
foreach var of varlist $tab_muni {
	
	local ytitle : var label `var'
				
		qui sum `var' if treat==1
		local mean1 : di  %4.2f r(mean)
		local sd1: di %4.2f r(sd)
		*local N1: di %4.0f r(N)
		
		qui sum `var' if treat==0
		local mean0 : di  %4.2f r(mean)
		local sd0: di %4.2f r(sd)  
		*local N0: di %4.0f r(N)
		
		ttest `var' , by(treat)
		local ttest: di %4.2f `r(p)'
		local diff: di %4.2f `mean1' - `mean0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
			
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\

	loc varN=`varN'+1		
}	
	*tex   \midrule
	*tex         N                 &   `N1' &    `N0'       &   \\  
	*tex   \midrule
texdoc close

keep if in_sample_reg_ecot == 1
keep if time==0

gen tipoe=1 if M1_103D==1
replace tipoe=0 if tipoe==. & M1_103D!=.
label var tipoe "Operator"

gen tipoe2=1 if M1_103D==2
replace tipoe2=0 if tipoe2==. & M1_103D!=.
label var tipoe2 "Supplier"

gen tipoe3=1 if M1_103D==3
replace tipoe3=0 if tipoe3==. & M1_103D!=.
label var tipoe3 "Beneficiary"

global main "tipoe tipoe2 tipoe3"
	
texdoc init "$outcomes/sum_stats_3b.tex" , replace force 
loc varN=1
cap mat drop coefs
foreach var of varlist $main {
	
	local ytitle : var label `var'

		qui sum `var' if treat==1 
		local mean1 : di  %4.2f r(mean)
		local sd1: di %4.2f r(sd)
		local N1: di %4.0f r(N)
		local N_total = r(N)
		
		qui sum `var' if treat==0 
		local mean0 : di  %4.2f r(mean)
		local sd0: di %4.2f r(sd)  
		local N0: di %4.0f r(N)
		local N_total: di %4.0f `N_total' + r(N)
		
		
		ttest `var' , by(treat)
		local ttest: di %4.2f `r(p)'
		local diff: di %4.2f `mean1' - `mean0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
			
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\
	
	if mod(`varN',3)==0 {
	*tex   \midrule
	*tex         N                 &   `N1' &    `N0'       &   `N_total'   \\  
	*tex   \midrule
	}
		
	loc varN=`varN'+1		
}	
texdoc close
restore
* Attrition: Tables/attrition_muni 



************ Table 2: Effects of eco-tourism promotion on number of tourists *** 
* output: Tables/ntourist_n
preserve

set obs 10000
global var "M1_412A ntourist M1_201A M1_201B"

* ANCOVA without controls 
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 if in_sample_ntourist == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var'_noc: reghdfe `var' treatXtime `var'_m1 if in_sample_ntourist == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var'_ANCOVA_noc: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA_nocont: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var'_ANCOVA: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT without controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat if time == 1 & in_sample_ntourist == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo RCT_`var'_noc: reghdfe `var' treatXtime if time == 1 & in_sample_ntourist == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var'_noc: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_nocont: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_ntourist == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo RCT_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_ntourist == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var': di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'
}
restore

* DiD without controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
gen randotreat_did = randotreat*time 
reghdfe `var' randotreat_did if in_sample_ntourist == 1, absorb(pair time) vce( cluster (codmpio)) 

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind randotreat_did
}

*I do the actual regression
eststo DiD_`var'_noc: reghdfe `var' treatXtime if in_sample_ntourist == 1, absorb(pair time) vce( cluster (codmpio)) 
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var'_noc: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_DiD_nocont: di %6.2f `r(N)'/`Nrandos'
}
restore

* DiD with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
gen randotreat_did = randotreat*time 
reghdfe `var' randotreat_did edu_anosXtime distanciaBogNMXtime DdistanciaBogm homicidiosXtime if in_sample_ntourist == 1, absorb(pair time) vce( cluster (codmpio)) 

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind randotreat_did
}

*I do the actual regression
eststo DiD_`var': reghdfe `var' treatXtime edu_anosXtime distanciaBogNMXtime DdistanciaBogm homicidiosXtime if in_sample_ntourist == 1, absorb(pair time) vce( cluster (codmpio)) 
qui sum `var' if treat==0 & time == 1 & in_sample_ntourist == 1
glo  mean_`var'_did: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_DiD: di %6.2f `r(N)'/`Nrandos'
}
restore

estout ancova_M1_412A_noc ancova_M1_412A ancova_ntourist_noc ancova_ntourist ancova_M1_201A_noc ancova_M1_201A ancova_M1_201B_noc ancova_M1_201B  using "$outcomes/ntourist_ANCOVA.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f ) labels ("Adjusted R$^2$")) keep(treatXtime) ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels(treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_M1_412A_ANCOVA_nocont} & ${rand_M1_412A_ANCOVA} & ${rand_ntourist_ANCOVA_nocont} & ${rand_ntourist_ANCOVA}  & ${rand_M1_201A_ANCOVA_nocont} & ${rand_M1_201A_ANCOVA} & ${rand_M1_201B_ANCOVA_nocont} & ${rand_M1_201B_ANCOVA} \\" ) 

estout  RCT_M1_412A_noc RCT_M1_412A RCT_ntourist_noc RCT_ntourist RCT_M1_201A_noc RCT_M1_201A RCT_M1_201B_noc RCT_M1_201B using "$outcomes/ntourist_RCT.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a  N  , fmt(%12.2f  %11.2gc ) labels ("Adjusted R$^2$" "Observations" )) keep(treatXtime) ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels(treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_M1_412A_RCT_nocont} & ${rand_M1_412A_RCT} & ${rand_ntourist_RCT_nocont} & ${rand_ntourist_RCT}  & ${rand_M1_201A_RCT_nocont} & ${rand_M1_201A_RCT} & ${rand_M1_201B_RCT_nocont} & ${rand_M1_201B_RCT} \\" ///
"Controls & No & Yes & No & Yes & No & Yes & No & Yes \\" ///
"Mean  Dep. Var. Control & ${mean_M1_412A_noc} & ${mean_M1_412A} & ${mean_ntourist_noc}  & ${mean_ntourist}  & ${mean_M1_201A_noc} & ${mean_M1_201A} & ${mean_M1_201B_noc} & ${mean_M1_201B} \\") 


estout  DiD_M1_412A_noc DiD_M1_412A DiD_ntourist_noc DiD_ntourist DiD_M1_201A_noc DiD_M1_201A DiD_M1_201B_noc DiD_M1_201B  using "$outcomes/ntourist_DiD.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a  N  , fmt(%12.2f  %11.2gc ) labels ("Adjusted R$^2$" "Observations" )) keep(treatXtime) ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels(treatXtime "After X Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_M1_412A_DiD_nocont} & ${rand_M1_412A_DiD} & ${rand_ntourist_DiD_nocont} & ${rand_ntourist_DiD} & ${rand_M1_201A_DiD_nocont} & ${rand_M1_201A_DiD} & ${rand_M1_201B_DiD_nocont} & ${rand_M1_201B_DiD}  \\" ///
"Controls & No & Yes & No & Yes & No & Yes & No & Yes \\" ///
"Mean  Dep. Var. Control  & ${mean_M1_412A_noc} & ${mean_M1_412A} & ${mean_ntourist_noc}  & ${mean_ntourist} & ${mean_M1_201A_noc} & ${mean_M1_201A} & ${mean_M1_201B_noc} & ${mean_M1_201B} \\") 
eststo clear


********************* Table 3: Spillovers using distance ***********************
* output: Tables/ntourist_spillovers

eststo spill_1: reghdfe ntourist treatXtime ntourist_m1 promotedshare28km edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_ntourist1: di %6.2f r(mean)

eststo spill_2: reghdfe ntourist treatXtime ntourist_m1 promotedshare28_56km edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_ntourist2: di %6.2f r(mean)

eststo spill_3: reghdfe ntourist treatXtime ntourist_m1 promotedshare56_133km edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_ntourist3: di %6.2f r(mean)

eststo spill_4: reghdfe ntourist treatXtime ntourist_m1 promotedshare28km promotedshare28_56km promotedshare56_133km edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_ntourist4: di %6.2f r(mean)

eststo spill_1a: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 promotedshare28km edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum LossSitio2Km if treat==0 & e(sample) == 1
glo  mean_defo1a: di %6.2f r(mean)

eststo spill_2a: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 promotedshare28_56km edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum LossSitio2Km if treat==0 & e(sample) == 1
glo  mean_defo2a: di %6.2f r(mean)

eststo spill_3a: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 promotedshare56_133km edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum LossSitio2Km if treat==0 & e(sample) == 1
glo  mean_defo3a: di %6.2f r(mean)

eststo spill_4a: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1  promotedshare28km promotedshare28_56km 	promotedshare56_133km  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS
qui sum LossCasa2Km if treat==0 & e(sample) == 1
glo  mean_defo4a: di %6.2f r(mean)

estout  spill_1 spill_2 spill_3 spill_4 spill_1a spill_2a spill_3a spill_4a using "$outcomes/ntourist_spillovers.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$")) keep(promotedshare28km promotedshare28_56km promotedshare56_133km treatXtime) ///
  cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels( treatXtime "Promoted" promotedshare28km "Promoted Share 28.1km" promotedshare28_56km "Promoted Share 28.1-56.5km" promotedshare56_133km "Promoted Share 56.5-133.1km") ///
  prefoot("\bottomrule" ///
"Mean  Dep. Var. Control  & ${mean_ntourist1}  & ${mean_ntourist2}  & ${mean_ntourist3}  & ${mean_ntourist4} &  ${mean_defo1a}  & ${mean_defo2a}  & ${mean_defo3a}  & ${mean_defo4b}\\") 
eststo clear

* ANCOVA *

************ Table 4: Effects of eco-tourism promotion on deforestation ********
preserve
set obs 10000
global var "LossSitio2Km AdefSitio2Km"

* ANCOVA with controls
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo anc_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
	qui sum `var'  if treat==0 & e(sample) & time == 1
	glo  `var'_mean: di %6.2f r(mean)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_ANC: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_ANC: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	di `p2'
	if `p2'< 1 {
	glo list_p_`var'_ANC: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_ANC: di %6.2f `p2_n'
	}	
}
restore

* ANCOVA with FOREST controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo anc_`var'_For : reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
	qui sum `var' if treat==0 & e(sample) & time == 1
	glo  `var'_mean_forest : di %6.2f r(mean)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_ANC_For: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_ANC_For: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_ANC_For: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_ANC_For: di %6.2f `p2_n'
	}
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_`var' == 1, absorb(pair) cluster(codmpio)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_RCT: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_RCT: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_RCT: di %6.2f `p2_n'
	}
}
restore

* RCT with FOREST controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if time == 1 & in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo rct_`var'_For: reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if time == 1 & in_sam_`var' == 1, absorb(pair) cluster(codmpio)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_RCT_For: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_RCT_For: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_RCT_For: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_RCT_For: di %6.2f `p2_n'
	}
}
restore

* ANCOVA with controls
preserve
set obs 10000
global var "LossCasa2Km AdefCasa2Km"
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	 qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_`var' == 1 & time == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo anc_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_`var' == 1 & time == 1, absorb(pair) cluster(codmpio)
	qui sum `var'  if treat==0 & e(sample) & time == 1
	glo  `var'_mean: di %6.2f r(mean)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_ANC: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_ANC: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_ANC: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_ANC: di %6.2f `p2_n'
	}
}
restore

* ANCOVA with FOREST controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1 & time == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo anc_`var'_For : reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1 & time == 1, absorb(pair) cluster(codmpio)
	qui sum `var' if treat==0 & e(sample) & time == 1
	glo  `var'_mean_forest : di %6.2f r(mean)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_ANC_For: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_ANC_For: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_ANC_For: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_ANC_For: di %6.2f `p2_n'
	}
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_`var' == 1, absorb(pair) cluster(codmpio)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_RCT: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_RCT: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_RCT: di %6.2f `p2_n'
	}
}
restore

* RCT with FOREST controls 
preserve
set obs 10000
foreach var of global var {
	capture drop randonum  terand randotreat diffbosque
	*Gen a variable where to store the treatment effects of the different randomization
	gen terand=.
	*Gen a variable where to store the treatment assignment of each randomization
	gen randotreat=.
	*Gen a variable to store the randon number that will guide the randomization
	gen randonum=.
	*Gen a variable to store the diff in forest to get similar randomizations
	gen diffbosque=.
	*decide how many randomizations
	local Nrandos=10000
	*set obs `Nrandos'
	*Start the loop of randomizations
	qui forval i =1/`Nrandos' {
		*Set the seed on each iteration to get different randomization
		set seed `i'
		*I gen an indicator by group (the strata of randomization) to indicate if they 
		*switch or maintain original treatment status
		*The conditional guarantees I generate one per Group
		sort pair
		replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
		*I copy thi for the rest of the group
		bys pair: egen randoind=max(randonum)
		*Mathematical trick to maintain or switch treatment status
		replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

		*Estimate the reg with the randomize treatment
		reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time == 1 & in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

		sum BosqueSitio2Km if randotreat==1 & e(sample)
		local foresttreat=`r(mean)'
		sum BosqueSitio2Km if randotreat==0 & e(sample)
		local forestcontrol=`r(mean)'

		*Store the ``treatment effect'' of the randomization
		replace terand=_b[randotreat]/_se[randotreat] in `i'
		replace diffbosque=`foresttreat'-`forestcontrol' in `i'
		drop randoind
	}

	*I do the actual regression
	eststo rct_`var'_For: reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time == 1 & in_sam_`var' == 1, absorb(pair) cluster(codmpio)

	*I count in how many cases I obtain a more adverse estimate
	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

	*Randomization inference p-value two-tails
	glo rand_`var'_RCT_For: di %6.2f `r(N)'/`Nrandos'

	qui sum diffbosque, d
	local bosquep25=`r(p25)'

	count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.)) & diffbosque<`bosquep25'
	local Nsimulbosquep25=`r(N)'

	count if  diffbosque<`bosquep25'
	local Nbosquep25=`r(N)'

	*Randomization inference p-value two-tails adjusting forest
	glo randadj_`var'_RCT_For: di %6.2f `Nsimulbosquep25'/`Nbosquep25'

	* Bonferroni hypothesis testing 
	local t = _b[treatXtime]/_se[treatXtime]
	local p =2*ttail(e(df_r),abs(`t'))
	local p2 =`p'*2
	if `p2'< 1 {
	glo list_p_`var'_RCT_For: di %6.2f `p2'
	} 
	else if `p2'>=1 {
	local p2_n = 1
	glo list_p_`var'_RCT_For: di %6.2f `p2_n'
	}
}
restore

**ANCOVA
estout anc_LossSitio2Km anc_LossSitio2Km_For anc_LossCasa2Km anc_LossCasa2Km_For anc_AdefSitio2Km anc_AdefSitio2Km_For anc_AdefCasa2Km anc_AdefCasa2Km_For using "$outcomes/forest_regs_backup.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f) labels ("Adjusted R$^2$" ))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"List (2019) p-value & ${list_p_LossSitio2Km_ANC}  & ${list_p_LossSitio2Km_ANC_For}  & ${list_p_LossCasa2Km_ANC}  & ${list_p_LossCasa2Km_ANC_For}  & ${list_p_AdefSitio2Km_ANC} & ${list_p_AdefSitio2Km_ANC_For} & ${list_p_AdefCasa2Km_ANC}  & ${list_p_AdefCasa2Km_ANC_For} \\" ///
"RI p-value & ${rand_LossSitio2Km_ANC}  & ${rand_LossSitio2Km_ANC_For}  & ${rand_LossCasa2Km_ANC}  & ${rand_LossCasa2Km_ANC_For}  & ${rand_AdefSitio2Km_ANC} & ${rand_AdefSitio2Km_ANC_For} & ${rand_AdefCasa2Km_ANC}  & ${rand_AdefCasa2Km_ANC_For}  \\" ///
"RI adj p-value & ${randadj_LossSitio2Km_ANC}  & ${randadj_LossSitio2Km_ANC_For}  & ${randadj_LossCasa2Km_ANC}  & ${randadj_LossCasa2Km_ANC_For}  & ${randadj_AdefSitio2Km_ANC} & ${randadj_AdefSitio2Km_ANC_For} & ${randadj_AdefCasa2Km_ANC}  & ${randadj_AdefCasa2Km_ANC_For}  \\" )

**RCT
estout rct_LossSitio2Km rct_LossSitio2Km_For rct_LossCasa2Km rct_LossCasa2Km_For rct_AdefSitio2Km rct_AdefSitio2Km_For rct_AdefCasa2Km rct_AdefCasa2Km_For  using "$outcomes/forest_regs_backup_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  r2_a N, fmt(%12.2f %11.2gc  ) labels ( "Adjusted R$^2$" "Observations" ))   ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"List (2019) p-value  & ${list_p_LossSitio2Km_RCT}  & ${list_p_LossSitio2Km_RCT_For}  & ${list_p_LossCasa2Km_RCT}  & ${list_p_LossCasa2Km_RCT_For}  & ${list_p_AdefSitio2Km_RCT} & ${list_p_AdefSitio2Km_RCT_For} & ${list_p_AdefCasa2Km_RCT}  & ${list_p_AdefCasa2Km_RCT_For} \\" ///
"RI p-value & ${rand_LossSitio2Km_RCT}  & ${rand_LossSitio2Km_RCT_For} & ${rand_LossCasa2Km_RCT}  & ${rand_LossCasa2Km_RCT_For}  & ${rand_AdefSitio2Km_RCT} & ${rand_AdefSitio2Km_RCT_For} & ${rand_AdefCasa2Km_RCT}  & ${rand_AdefCasa2Km_RCT_For}  \\" ///
"RIadj p-value & ${randadj_LossSitio2Km_RCT}  & ${randadj_LossSitio2Km_RCT_For} & ${randadj_LossCasa2Km_RCT}  & ${randadj_LossCasa2Km_RCT_For}  & ${randadj_AdefSitio2Km_RCT} & ${randadj_AdefSitio2Km_RCT_For} & ${randadj_AdefCasa2Km_RCT}  & ${randadj_AdefCasa2Km_RCT_For}  \\" ///
"Forest Control & No & Yes & No & Yes & No & Yes & No & Yes \\" ///
"Radius & 2 & 2 & 2 & 2 & 2 & 2 & 2 & 2 \\" ///
"Mean  Dep. Var. Control  & ${LossSitio2Km_mean}  & ${LossSitio2Km_mean}  & ${LossCasa2Km_mean}  & ${LossCasa2Km_mean}   & ${AdefSitio2Km_mean}  & ${AdefSitio2Km_mean} & ${AdefCasa2Km_mean} & ${AdefCasa2Km_mean} \\"  ) 
eststo clear


******* Table 5: Effects of eco-tourism promotion on businesses' outcomes ******
* Variables nworkers M1_205

preserve
set obs 10000
global var "nworkers M1_205 M1_312" 
* ANCOVA without controls 
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 if in_sample_busoutcomes == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var'_nocont: reghdfe `var' treatXtime `var'_m1 if in_sample_busoutcomes == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA_nocont: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_ANCOVA_nocont: di %6.2f `p2'
}
restore

* ANCOVA with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_busoutcomes == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_busoutcomes == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_ANCOVA: di %6.2f `p2'
}
restore

* RCT without controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat if time == 1 & in_sample_busoutcomes == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_nocont: reghdfe `var' treatXtime if time == 1 & in_sample_busoutcomes == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_nocont: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_RCT_nocont: di %6.2f `p2'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_RCT: di %6.2f `p2'
}
restore

preserve
set obs 10000
global var "pobreza_inc tourism_income month_income_all agro_income"

* ANCOVA without controls 
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 if in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var'_nc: reghdfe `var' treatXtime `var'_m1 if in_sample_reg_ecot == 1, absorb(pair) cluster(codmpio)
qui sum `var'  if treat==0 & e(sample) & time == 1
glo  `var' _mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA_nc: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_ANC_nc: di %6.2f `p2'
}
restore

* ANCOVA with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_ANC: di %6.2f `p2'
}
restore

* RCT without controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat if time == 1 & in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_nc: reghdfe `var' treatXtime if time == 1 & in_sample_reg_ecot == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_nc: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_RCT_nc: di %6.2f `p2'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_reg_ecot == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'

* Bonferroni hypothesis testing 
local t = _b[treatXtime]/_se[treatXtime]
local p =2*ttail(e(df_r),abs(`t'))
local p2 =`p'*2
glo list_p_`var'_RCT: di %6.2f `p2'
}
restore

* ANCOVA 
* output: Tables/others_outcomes  
estout  ancova_nworkers ancova_M1_205 ancova_M1_312 ancova_pobreza_inc ancova_tourism_income ancova_agro_income ancova_month_income_all using "$outcomes/others_outcomes.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f ) labels  ("Adjusted R$^2$") )  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") prefoot("\bottomrule" ///
"List(2019) p-value & ${list_p_nworkers_ANCOVA} & ${list_p_M1_205_ANCOVA} & ${list_p_M1_312_ANCOVA} & ${list_p_pobreza_inc_ANC} & ${list_p_tourism_income_ANC} & ${list_p_agro_income_ANC} & ${list_p_month_income_all_ANC} \\" ///
"RI p-value & ${rand_nworkers_ANCOVA} & ${rand_M1_205_ANCOVA} & ${rand_M1_312_ANCOVA} & ${rand_pobreza_inc_ANCOVA} & ${rand_tourism_income_ANCOVA} & ${rand_agro_income_ANCOVA} & ${rand_month_income_all_ANCOVA} \\" )

* output: Tables/others_outcomes_nocont
estout ancova_nworkers_nocont ancova_M1_205_nocont ancova_M1_312_nocont ancova_pobreza_inc_nc ancova_tourism_income_nc ancova_agro_income_nc ancova_month_income_all_nc using "$outcomes/others_outcomes_nocont.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f ) labels  ("Adjusted R$^2$") )  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") prefoot("\bottomrule" ///
"List(2019) p-value & ${list_p__nworkers_ANCOVA_nocont} & ${list_p_M1_205_ANCOVA_nocont} & ${list_p_M1_312_ANCOVA_nocont} & ${list_p_pobreza_inc_ANCOVA_nc} & ${list_p_tourism_income_ANC_nc} & ${list_p_agro_income_ANC_nc} & ${list_p_month_income_all_ANC_nc} \\" ///
"RI p-value & ${rand_nworkers_ANCOVA_nocont} & ${rand_M1_205_ANCOVA_nocont} & ${rand_M1_312_ANCOVA_nocont} & ${rand_pobreza_inc_ANCOVA_nc} & ${rand_tourism_income_ANCOVA_nc} & ${rand_agro_income_ANCOVA_nc}  & ${rand_month_income_all_ANCOVA_nc} \\" )

* RCT:  
* output: Tables/others_outcomes_rct 
estout rct_nworkers rct_M1_205 rct_M1_312 rct_pobreza_inc rct_tourism_income rct_agro_income rct_month_income_all using "$outcomes/others_outcomes_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(N r2_a , fmt(%11.2gc %12.2f  ) labels ("Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"List(2019) p-value & ${list_p_nworkers_RCT} & ${list_p_M1_205_RCT} & ${list_p_M1_312_RCT} & ${list_p_pobreza_inc_RCT} & ${list_p_tourism_income_RCT} & ${list_p_agro_income_RCT} & ${list_p_month_income_all_RCT} \\" ///
"RI p-value & ${rand_nworkers_RCT} & ${rand_M1_205_RCT}  & ${rand_M1_312_RCT} & ${rand_pobreza_inc_RCT} & ${rand_tourism_income_RCT} & ${rand_agro_income_RCT}  & ${rand_month_income_all_RCT}  \\" ///
"Controls & Yes & Yes & Yes & Yes & Yes & Yes & Yes \\" ///
"Mean  Dep. Var. Control  & ${nworkers_mean }  & ${M1_205_mean}  & ${M1_312_mean} & ${pobreza_inc_mean }  & ${tourism_income_mean } & ${agro_income_mean } & ${month_income_all_mean}  \\") 

* output: Tables/others_outcomes_nocont_rct 
estout  rct_nworkers_nocont rct_M1_205_nocont rct_M1_312_nocont rct_pobreza_inc_nc rct_tourism_income_nc rct_agro_income_nc rct_month_income_all_nc using "$outcomes/others_outcomes_nocont_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(N r2_a , fmt(%11.2gc %12.2f  ) labels ("Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"List(2019) p-value & ${list_p_nworkers_RCT_nocont} & ${list_p_M1_205_RCT_nocont} & ${list_p_M1_312_RCT_nocont} & ${list_p_pobreza_inc_RCT_nc} & ${list_p_tourism_income_RCT_nc} & ${list_p_agro_income_RCT_nc} & ${list_p_month_income_all_RCT_nc} \\" ///
"RI p-value & ${rand_nworkers_RCT_nocont} & ${rand_M1_205_RCT_nocont} & ${rand_M1_312_RCT_nocont} &  ${rand_pobreza_inc_RCT_nc} & ${rand_tourism_income_RCT_nc} & ${rand_agro_income_RCT_nc}   & ${rand_month_income_all_RCT_nc}\\" ///
"Controls & No & No & No & No & No & No & No \\" ///
"Mean  Dep. Var. Control  & ${nworkers_mean } & ${M1_205_mean} & ${M1_312_mean} & ${pobreza_inc_mean }  & ${tourism_income_mean } & ${agro_income_mean } & ${month_income_all_mean}  \\") 
eststo clear 


********************************************************************************
******************************* APPENDIX ***************************************
********************************************************************************

************** Table A.1: Effects of eco-tourism promotion on number of tourists using Diff-in-Diff specification
* Computed in T2

************ Table A.2: Heterogeneous results **********************************
* output: Tables/ntourist_het
eststo est1:  reghdfe ntourist treatXtime ntourist_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   cluster(codmpio) //ANCOVA CONTROLS
gen sample_heterog = 1 if e(sample)
qui sum ntourist if treat==0 & sample_heterog==1 
glo tourist1: di %6.2f r(mean)
drop sample_heterog

eststo est2: reghdfe ntourist  treatX_Exist treatX_NewProm treatX_NewCha ntourist_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)  cluster(codmpio)
gen sample_heterog = 1 if e(sample)
qui sum ntourist if treat==0 & sample_heterog==1 
glo tourist2: di %6.2f r(mean)
drop sample_heterog

eststo est3: reghdfe ntourist  viox* ntourist_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)  cluster(codmpio)
gen sample_heterog = 1 if e(sample)
qui sum ntourist if treat==0 & sample_heterog==1 
glo tourist3: di %6.2f r(mean)
drop sample_heterog

eststo est4: reghdfe M1_412A treatXtime M1_412A_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if M1_103D!=3, absorb(pair)   cluster(codmpio) //ANCOVA CONTROLS
gen sample_heterog = 1 if e(sample)
qui sum M1_412A if treat==0 & sample_heterog==1 
glo M1_412A_1: di %6.2f r(mean)

eststo est5: reghdfe M1_412A  treatX_Exist treatX_NewProm treatX_NewCha M1_412A_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if sample_heterog == 1, absorb(pair)  cluster(codmpio)
qui sum M1_412A if treat==0 & sample_heterog==1 
glo M1_412A_2: di %6.2f r(mean)

eststo est6: reghdfe M1_412A  viox* M1_412A_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if sample_heterog == 1, absorb(pair)  cluster(codmpio)
qui sum M1_412A if treat==0 & sample_heterog==1 
glo M1_412A_3: di %6.2f r(mean)
drop sample_heterog

estout est1 est2 est3 est4 est5 est6  using "${outcomes}/ntourist_het.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime treatX* viox*) varlabels( treatXtime "Promoted" treatX_Exist "Promoted X Existing" treatX_NewProm "Promoted X Promising"  treatX_NewCha "Promoted X Challenging" viox_bajo "Promoted X Low" viox_medio "Promoted X Medium" viox_alto "Promoted X High") ///
prefoot("\bottomrule" ///
"Mean  Dep. Var. Control  & ${tourist1}  & ${tourist2}  & ${tourist3}    & ${M1_412A_1} & ${M1_412A_2} & ${M1_412A_3}  \\") 
eststo clear

************************** Table A.3: Robustness of spillover results *************************
 ** Ntourists spillovers using Heb **

eststo spill_1a_tourist: reghdfe ntourist treatXtime ntourist_m1 treatedngb28km_area munisample28km_area ///
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers until 28Km
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_defo1a_tourist: di %6.2f r(mean)

eststo spill_2a_tourist: reghdfe ntourist treatXtime ntourist_m1 treatngb28_56km_area munisamp28_56km_area /// 
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers from 28 to 56Km 
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_defo2a_tourist: di %6.2f r(mean)

eststo spill_3a_tourist: reghdfe ntourist treatXtime ntourist_m1 treatngb56_133km_area munisamp56_133km_area /// 
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers from 56 to 133Km 
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_defo3a_tourist: di %6.2f r(mean)

eststo spill_4a_tourist: reghdfe ntourist treatXtime ntourist_m1 treatedngb28km munisample28km_area ///
                         treatngb28_56km_area munisamp28_56km_area  ///
						 treatngb56_133km_area munisamp56_133km_area  ///
						 edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_ntourist == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + All previous spillovers
qui sum ntourist if treat==0 & in_sample_ntourist == 1 & time == 1
glo  mean_defo4a_tourist: di %6.2f r(mean)

 ** Deforestation spillovers using Heb ** 
eststo spill_1a_sitio: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 treatedngb28km_area munisample28km_area ///
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers until 28Km
qui sum LossSitio2Km if treat==0 & e(sample) == 1 
glo  mean_defo1a_sitio: di %6.2f r(mean)

eststo spill_2a_sitio: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 treatngb28_56km_area munisamp28_56km_area /// 
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers from 28 to 56Km 
qui sum LossSitio2Km if treat==0 & e(sample) == 1 
glo  mean_defo2a_sitio: di %6.2f r(mean)

eststo spill_3a_sitio: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 treatngb56_133km_area munisamp56_133km_area /// 
                  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + Spillovers from 56 to 133Km 
qui sum LossSitio2Km if treat==0 & e(sample) == 1 
glo  mean_defo3a_sitio: di %6.2f r(mean)

eststo spill_4a_sitio: reghdfe LossSitio2Km treatXtime LossSitio2Km_m1 treatedngb28km munisample28km_area ///
                         treatngb28_56km_area munisamp28_56km_area  ///
						 treatngb56_133km_area munisamp56_133km_area ///
						 edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair)   vce(robust) //ANCOVA CONTROLS + All previous spillovers
qui sum LossSitio2Km  if treat==0 & e(sample) == 1 
glo  mean_defo4a_sitio: di %6.2f r(mean)

estout  spill_1a_tourist spill_2a_tourist spill_3a_tourist spill_4a_tourist spill_1a_sitio spill_2a_sitio spill_3a_sitio spill_4a_sitio using "$outcomes/deforestation&ntourist_heb2021_spillovers.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$")) keep(treatXtime treatedngb28km_area treatngb28_56km_area treatngb56_133km_area munisample28km_area munisamp28_56km_area munisamp56_133km_area) ///
  cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels( treatXtime "Promoted" treatedngb28km_area "\% Promoted Area 28.1km" treatngb28_56km_area "\% Promoted Area 28.1-56.5km" treatngb56_133km_area "\% Promoted Area 56.5-133.1km" munisample28km_area "\% Area of munis in sample 28.1km" munisamp28_56km_area "\% Area of munis in sample 28.1-56.5km" munisamp56_133km_area "\% Area of munis in sample 56.5-133.1km" ) ///
  prefoot("\bottomrule" ///
"Mean  Dep. Var. Control  & ${mean_defo1a_tourist} & ${mean_defo2a_tourist} & ${mean_defo3a_tourist} & ${mean_defo4a_tourist} & ${mean_defo1a_sitio}  & ${mean_defo2a_sitio}  & ${mean_defo3a_sitio}   & ${mean_defo4a_sitio} \\") 
eststo clear

*********** Table A.4: Effects of eco-tourism promotion on socio-economic variables without extra surveys
************************* without extra surveys ********************************
* output: Tables/others_outcomes_nw
* output: Tables/others_outcomes_rct_nw

global var "M1_412A"
foreach var of global var {
	
preserve
eststo est_`var':  reghdfe `var'  treatXtime `var'_m1   edu_anos distanciaBogNM  DdistanciaBogm homicidios if M1_103D!=3 & reg_m==1 , absorb(pair)   cluster(codmpio) //ANCOVA CONTROLS
gen sample_= 1 if e(sample)
qui sum `var' if treat == 0 & sample_ == 1 
glo `var': di %6.2f r(mean)

keep if time==1 
eststo est_`var'_rct: reghdfe `var'  treatXtime  edu_anos distanciaBogNM  DdistanciaBogm homicidios if M1_103D!=3 & sample_==1 & reg_m==1 , absorb(pair) cluster(codmpio) ///RCT CONTROLS

restore
}

global var " ntourist nworkers M1_205 pobreza_inc tourism_income month_income_all"
foreach var of global var {

preserve
eststo est_`var':  reghdfe `var'  treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm homicidios if  reg_m==1   , absorb(pair)   cluster(codmpio)  //ANCOVA CONTROLS 
gen sample_ = 1 if e(sample)
qui sum `var' if treat == 0 & e(sample) == 1
glo `var': di %6.2f r(mean)


keep if time == 1 
eststo est_`var'_rct: reghdfe `var'  treatXtime  edu_anos distanciaBogNM  DdistanciaBogm homicidios if sample_==1 & reg_m==1 , absorb(pair) cluster(codmpio)  //RCT CONTROLS
qui sum `var' if treat == 0 & e(sample) == 1
glo `var'_rct: di %6.2f r(mean)
drop sample_

restore
}

* ANCOVA 
* output: Tables/others_outcomes_nw
estout est_ntourist est_M1_412A est_nworkers est_M1_205 est_pobreza_inc est_tourism_income est_month_income_all  using "${outcomes}/others_outcomes_nw.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f) labels ("Adjusted R$^2$" ))   ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels(treatXtime "Promoted") prefoot("\bottomrule")

***RCT
* output: Tables/others_outcomes_rct_nw
estout est_ntourist_rct est_M1_412A_rct est_nworkers_rct est_M1_205_rct est_pobreza_inc_rct est_tourism_income_rct est_month_income_all_rct using "${outcomes}/others_outcomes_rct_nw.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels(treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"Mean  Dep. Var. Control & ${ntourist_rct}  & ${M1_412A}  & ${nworkers_rct}  & ${M1_205_rct} & ${pobreza_inc_rct } & ${tourism_income_rct}  & ${month_income_all_rct} \\") 
eststo clear

*********** Table A.5: Heterogeneous results for others outcomes ***************
* output: Tables/ntourist_het_nw
global var "bosqprop2Km"
foreach var of global var {

eststo `var':  reghdfe `var'  treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_`var'== 1, absorb(pair) cluster(codmpio) //ANCOVA CONTROLS
gen sample_=1 if e(sample)
qui sum `var' if treatXtime==0 & sample_==1 
glo `var'_1: di %6.2f r(mean)
drop sample_

eststo `var'_interact: reghdfe `var'  treatX_Exist treatX_NewProm treatX_NewCha D_Exit D_NewProm D_NewCha `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_`var'== 1, absorb(pair) cluster(codmpio) //ANCOVA CONTROLS + INTERACT
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'_2: di %6.2f r(mean)
drop sample_
}

global var "bosqprop2KmC"
foreach var of global var {

eststo `var':  reghdfe `var'  treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_`var'== 1, absorb(pair) cluster(codmpio) //ANCOVA CONTROLS
gen sample_=1 if e(sample)
qui sum `var' if treatXtime==0 & sample_==1 
glo `var'_1: di %6.2f r(mean)
drop sample_

eststo `var'_interact2: reghdfe `var' typex_operador typex_proveedor typex_beneficiario D_operator D_supplier D_beneficiary `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_`var'== 1, absorb(pair) cluster(codmpio) //ANCOVA CONTROLS + INTERACT
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'_2: di %6.2f r(mean)
drop sample_
}

global var "nworkers"

foreach var of global var {

eststo `var':  reghdfe `var'  treatXtime `var'_m1   edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio) //ANCOVA CONTROLS
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'1: di %6.2f r(mean)
drop sample_

eststo `var'_interact: reghdfe `var'  treatX_Exist treatX_NewProm treatX_NewCha D_Exit D_NewProm D_NewCha `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)  cluster(codmpio)
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'2: di %6.2f r(mean)
drop sample_

eststo `var'_interact2: reghdfe `var'  typex_operador typex_proveedor D_operator D_supplier `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)  cluster(codmpio)
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'3: di %6.2f r(mean)
drop sample_
}

global var "month_income_all"

foreach var of global var {

eststo `var':  reghdfe `var'  treatXtime `var'_m1   edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)   cluster(codmpio) //ANCOVA CONTROLS
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'1: di %6.2f r(mean)
drop sample_

eststo `var'_interact: reghdfe `var'  treatX_Exist treatX_NewProm treatX_NewCha D_Exit D_NewProm D_NewCha `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)  cluster(codmpio)
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'2: di %6.2f r(mean)
drop sample_

eststo `var'_interact2: reghdfe `var' typex_operador typex_proveedor typex_beneficiario D_operator D_supplier D_beneficiary `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios if in_sample_reg_ecot == 1, absorb(pair)  cluster(codmpio)
gen sample_=1 if e(sample)
qui sum `var' if treat==0 & sample_==1 
glo `var'3: di %6.2f r(mean)
drop sample_
}

estout bosqprop2Km_interact nworkers_interact month_income_all_interact using "${outcomes}/ntourist_het_nw_panelA.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatX*) varlabels(treatX_Exist "Promoted X Existing" treatX_NewProm "Promoted X Promising"  treatX_NewCha "Promoted X Challenging" typex_operador  "Promoted X Operator" typex_proveedor "Promoted X Supplier" type2x_beneficiariod "Promoted X Beneficiary Direct" type2x_beneficiarioin "Promoted X Beneficiary Indirect" typex_beneficiario  "Promoted X Beneficiary ") ///
prefoot("\bottomrule" ///
"Mean  Dep. Var. Control  & ${bosqprop2Km_2}  & ${nworkers2}  & ${month_income_all2} \\")

estout bosqprop2KmC_interact2 nworkers_interact2 month_income_all_interact2 using "${outcomes}/ntourist_het_nw_panelB.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  N r2_a , fmt(%11.2gc %12.2f ) labels ( "Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(typex_*) varlabels(treatX_Exist "Promoted X Existing" treatX_NewProm "Promoted X Promising"  treatX_NewCha "Promoted X Challenging" typex_operador  "Promoted X Operator" typex_proveedor "Promoted X Supplier" type2x_beneficiariod "Promoted X Beneficiary Direct" type2x_beneficiarioin "Promoted X Beneficiary Indirect" typex_beneficiario  "Promoted X Beneficiary ") ///
prefoot("\bottomrule" ///
"Mean  Dep. Var. Control  & ${bosqprop2KmC_2}  & ${nworkers3} & ${month_income_all3} \\")

eststo clear

************ Table A.6: Robustness to the dependent variable measurement **********
*** ANCOVA regressions ***
global var "LossSitio2Km LossCasa2Km"
foreach var of global var {
eststo anc_dummy_`var': reghdfe D_`var' D_`var'_m1 treatXtime   edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)

eststo anc_hs_`var': reghdfe hs_`var' hs_`var'_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)

eststo anc_ln_`var': reghdfe ln_`var' ln_`var'_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)
}

** Simulations **
* ANCOVA with forest controls Dummy *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe D_`var' randotreat D_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_D_`var': reghdfe D_`var' treatXtime D_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum D_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_D_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_D: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with forest controls HS *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe hs_`var' randotreat hs_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_hs_`var': reghdfe hs_`var' treatXtime hs_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum hs_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_hs_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_hs: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with forest controls logaritmic *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe ln_`var' randotreat ln_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_ln_`var': reghdfe ln_`var' treatXtime ln_`var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum ln_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_ln_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_ln: di %6.2f `r(N)'/`Nrandos'
}
restore
**ANCOVA
estout anc_D_LossSitio2Km anc_hs_LossSitio2Km anc_ln_LossSitio2Km anc_D_LossCasa2Km anc_hs_LossCasa2Km anc_ln_LossCasa2Km using  "$outcomes/forest_regs_hyperbolicsine_anc.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f) labels ("Adjusted R$^2$" ))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_LossSitio2Km_ANC_D}  & ${rand_LossSitio2Km_ANC_hs}  & ${rand_LossSitio2Km_ANC_ln}  & ${rand_LossCasa2Km_ANC_D} & ${rand_LossCasa2Km_ANC_hs} & ${rand_LossCasa2Km_ANC_ln} \\" )



*** RCT regressions ***
eststo rct_dummy_1: reghdfe D_LossSitio2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)

eststo rct_hs_1: reghdfe hs_LossSitio2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)

eststo rct_ln_1: reghdfe ln_LossSitio2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)

eststo rct_dummy_2: reghdfe D_LossCasa2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossCasa2Km == 1, absorb(pair) cluster(codmpio)

eststo rct_hs_2: reghdfe hs_LossCasa2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossCasa2Km == 1, absorb(pair) cluster(codmpio)

eststo rct_ln_2: reghdfe ln_LossCasa2Km treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time ==1 & in_sam_LossCasa2Km == 1, absorb(pair) cluster(codmpio)

** Simulations **
* RCT with forest controls Dummy *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe D_`var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo RCT_D_`var': reghdfe D_`var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum D_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_D_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_D: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with forest controls HS *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe hs_`var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo RCT_hs_`var': reghdfe hs_`var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum hs_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_hs_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_hs: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with forest controls logaritmic *
preserve
set obs 10000
* This code is using hiperbolice sine variable 
global var "LossSitio2Km LossCasa2Km"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe ln_`var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
} 

*I do the actual regression
eststo RCT_ln_`var': reghdfe ln_`var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sam_`var' == 1, absorb(pair) cluster(codmpio)
qui sum ln_`var'  if treat==0 & e(sample) & time == 1
glo  `var'_ln_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_ln: di %6.2f `r(N)'/`Nrandos'
}
restore

estout RCT_D_LossSitio2Km RCT_hs_LossSitio2Km RCT_ln_LossSitio2Km RCT_D_LossCasa2Km RCT_hs_LossCasa2Km RCT_ln_LossCasa2Km using "$outcomes/forest_regs_hyperbolicsine_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  r2_a N, fmt(%12.2f %11.2gc  ) labels ( "Adjusted R$^2$" "Observations" ))   ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_LossSitio2Km_RCT_D}  & ${rand_LossSitio2Km_RCT_hs} & ${rand_LossSitio2Km_RCT_ln}  & ${rand_LossCasa2Km_RCT_D}  & ${rand_LossCasa2Km_RCT_hs} & ${rand_LossCasa2Km_RCT_ln}\\" ///
"Forest Control & Yes  & Yes & Yes & Yes & Yes \\" ///
"Radius & 2 & 2 & 2 & 2 & 2 & 2 \\" ///
"Mean  Dep. Var. Control  & ${LossSitio2Km_D_mean}  & ${LossSitio2Km_hs_mean}  & ${LossSitio2Km_ln_mean}  & ${LossCasa2Km_D_mean}   & ${LossCasa2Km_hs_mean}  & ${LossCasa2Km_ln_mean} \\")
eststo clear

************* Table A.7: Effects of eco-tourism promotion on number of workers *********************************

* ANCOVA
eststo ancova_nwor1: reghdfe nworkers nworkers_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<12, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_ancova_nwor1: di %6.2f r(mean)
eststo ancova_nwor2: reghdfe nworkers nworkers_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<6, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_ancova_nwor2: di %6.2f r(mean)
eststo ancova_nwor3: reghdfe nworkers nworkers_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<4, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_ancova_nwor3: di %6.2f r(mean)

*RCT
eststo rct_nwor1: reghdfe nworkers treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<12, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_rct_nwor1: di %6.2f r(mean)
eststo rct_nwor2: reghdfe nworkers treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<6, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_rct_nwor2: di %6.2f r(mean)
eststo rct_nwor3: reghdfe nworkers treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_busoutcomes == 1 & nworkers<4, absorb(pair) cluster(codmpio)
qui sum nworkers  if treat==0 & e(sample) & time == 1
glo  m_rct_nwor3: di %6.2f r(mean)

* output: Tables/nworkers_robust 
estout ancova_nwor1 ancova_nwor2 ancova_nwor3 using "$outcomes/nworkers_robust.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f ) labels ( "Adjusted R$^2$") )  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") prefoot("\bottomrule" )

* RCT:  
* output: Tables/others_outcomes_rct
estout rct_nwor1 rct_nwor2 rct_nwor3 using "$outcomes/nworkers_robust_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats( N r2_a , fmt(%11.2gc %12.2f  ) labels ("Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"Controls & Yes & Yes & Yes \\" ///
"Mean  Dep. Var. Control  & ${m_ancova_nwor1 }  & ${m_ancova_nwor2 }  & ${m_ancova_nwor3 } \\")	
eststo clear 

********************************************************************************
******************************* APPENDIX ONLINE ********************************
********************************************************************************

************ Table B.1: Summary attractions ************************************
* output: Tables/atractions_eng
preserve
keep if unique_id_sitio == 1 & time == 0
tab M1_109, sort  // Check order are the same
estpost tabulate M1_109_eng if M1_109_eng != ".", sort 
*esttab using "$outcomes/atractions_eng.tex", cells("pct(f(2))") noobs replace
restore

************ Table B.2: Paired municipalites ***********************************
preserve
merge m:1 codmpio using "$base_out/Temporary/distance_codmpio.dta"
keep M1_103A M1_103B codmpio pair Direct_distKm treat
duplicates drop codmpio, force
merge 1:1 codmpio using "$base_out/Temporary/Grupo_ordered.dta"
sort ordered_group depto name_mpio

keep ordered_group depto name_mpio codmpio Direct_distKm

order(depto name_mpio ordered_group codmpio Direct_distKm)
replace Direct_distKm=round(Direct_distKm, 0.10)
texsave using "$base_out/Temporary/tabla_parejas_v2.tex",  replace

restore
************ Table B.3: Summary statistics forest cover ************************
preserve
reghdfe LossSitio2Km LossSitio2Km_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km, absorb(pair) cluster(codmpio)
gen in_sample_Sitio2Km = 1 if e(sample) 
bys NIM: egen in_sample_forest_siti = max(in_sample_Sitio2Km)

reghdfe LossCasa2Km LossCasa2Km_m1 treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km, absorb(pair) cluster(codmpio)
gen in_sample_Casa2Km = 1 if e(sample)
bys NIM: egen in_sample_forest = max(in_sample_Casa2Km)

global table1 "BosqueCasa1Km BosqueCasa2Km BosqueCasa5Km"
global table2 "BosqueSitio1Km BosqueSitio2Km BosqueSitio5Km"

keep if time == 0

texdoc init "$outcomes/sum_stats_forest.tex" , replace force 
loc varN=1
cap mat drop coefs
foreach var of varlist $table1 {
	
	local ytitle : var label `var'
				
		qui sum `var' if treat==1 & in_sample_forest == 1
		local mean1 : di  %15.2fc r(mean)
		local m1 = r(mean)
		local sd1: di %15.2fc r(sd)
		local N1: di %15.0fc r(N)
		
		qui sum `var' if treat==0 & in_sample_forest == 1
		local mean0 : di  %15.2fc r(mean)
		local m0 = r(mean)
		local sd0: di %15.2fc r(sd)  
		local N0: di %15.0f r(N)
		
		ttest `var' if in_sample_forest == 1, by(treat)
		local ttest: di %15.2f `r(p)'
		local diff: di %15.2fc `m1' - `m0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
	
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\
	if mod(`varN',3)==0 {
	tex   \midrule
	tex         N                 &   `N1' &    `N0'       &         \\  
	tex   \midrule
	}
	loc varN=`varN'+1	
}	

foreach var of varlist $table2 {
	
	local ytitle : var label `var'
				
		qui sum `var' if treat==1 & in_sample_forest_siti == 1
		local mean1 : di  %15.2fc r(mean)
		local m1 = r(mean)
		local sd1: di %15.2fc r(sd)
		local N1: di %15.0fc r(N)
		
		qui sum `var' if treat==0 & in_sample_forest_siti == 1
		local mean0 : di  %15.2fc r(mean)
		local m0 = r(mean)
		local sd0: di %15.2fc r(sd)  
		local N0: di %15.0f r(N)
		
		ttest `var' if in_sample_forest_siti == 1, by(treat)
		local ttest: di %15.2f `r(p)'
		local diff: di %15.2fc `m1' - `m0'
		loc star= ""
		if ((`r(p)' < 0.1) )  loc star = "*" 
		if ((`r(p)' < 0.05) ) loc star = "**" 
		if ((`r(p)' < 0.01) ) loc star = "***" 	
	
	tex \parbox[l]{7cm}{`ytitle'} & `mean1' & `mean0' & `diff'`star'\\
	tex                          & (`sd1') & (`sd0' ) & [`ttest'] \\
	if mod(`varN',3)==0 {
	tex   \midrule
	tex         N                 &   `N1' &    `N0'       &         \\  
	tex   \midrule
	}
	loc varN=`varN'+1	
}

texdoc close
restore

************* Table B.4: Effects of eco-tourism promotion on deforestation rate *****

preserve
set obs 10000
global var "LossSit5Km_r"

* ANCOVA with controls
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm  homicidios if in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm  homicidios if in_sample_`var' == 1, absorb(pair) cluster(codmpio)
qui sum `var'  if treat==0 & e(sample) & time == 1
glo  `var'_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm  homicidios BosqueSitio5Km if in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var'_Forest: reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm  homicidios BosqueSitio5Km if in_sample_`var' == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM  DdistanciaBogm  homicidios if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM  DdistanciaBogm  homicidios if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM  DdistanciaBogm  homicidios BosqueSitio5Km if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_Forest: reghdfe `var' treatXtime edu_anos distanciaBogNM  DdistanciaBogm  homicidios BosqueSitio5Km if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

preserve
set obs 10000
global var "LossSit2Km_r bosqprop2Km"

* ANCOVA with controls
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_`var' == 1, absorb(pair) cluster(codmpio)
qui sum `var'  if treat==0 & e(sample) & time == 1
glo  `var'_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var'_Forest : reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if in_sample_`var' == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean_forest : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_Forest: reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueSitio2Km if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

preserve
set obs 10000
global var "bosqprop2KmC"

* ANCOVA with controls
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_`var' == 1 & time == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sample_`var' == 1 & time == 1, absorb(pair) cluster(codmpio)
qui sum `var'  if treat==0 & e(sample) & time == 1
glo  `var'_mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sample_`var' == 1 & time == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo anc_`var'_Forest : reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if in_sample_`var' == 1 & time == 1, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean_forest : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANC_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with FOREST controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time == 1 & in_sample_`var' == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_Forest: reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios BosqueCasa2Km if time == 1 & in_sample_`var' == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_Forest: di %6.2f `r(N)'/`Nrandos'
}
restore

**ANCOVA
estout anc_LossSit2Km_r anc_LossSit2Km_r_Forest anc_LossSit5Km_r anc_LossSit5Km_r_Forest anc_bosqprop2Km anc_bosqprop2Km_Forest anc_bosqprop2KmC anc_bosqprop2KmC_Forest using "$outcomes/forest_regs.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f) labels ("Adjusted R$^2$" ))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_LossSit2Km_r_ANC}  & ${rand_LossSit2Km_r_ANC_Forest}  & ${rand_LossSit5Km_r_ANC}  & ${rand_LossSit5Km_r_ANC_Forest}  & ${rand_bosqprop2Km_ANC} & ${rand_bosqprop2Km_ANC_Forest} & ${rand_bosqprop2KmC_ANC}  & ${rand_bosqprop2KmC_ANC_Forest}  \\" )

**RCT
estout rct_LossSit2Km_r rct_LossSit2Km_r_Forest rct_LossSit5Km_r rct_LossSit5Km_r_Forest rct_bosqprop2Km rct_bosqprop2Km_Forest rct_bosqprop2KmC rct_bosqprop2KmC_Forest  using "$outcomes/forest_regs_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(  r2_a N, fmt(%12.2f %11.2gc  ) labels ( "Adjusted R$^2$" "Observations" ))   ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_LossSit2Km_r_RCT}  & ${rand_LossSit2Km_r_RCT_Forest} & ${rand_LossSit5Km_r_RCT}  & ${rand_LossSit5Km_r_RCT_Forest}  & ${rand_bosqprop2Km_RCT} & ${rand_bosqprop2Km_RCT_Forest} & ${rand_bosqprop2KmC_RCT}  & ${rand_bosqprop2KmC_RCT_Forest}  \\" ///
"Forest Control & No & Yes & No & Yes & No & Yes & No & Yes \\" ///
"Radius & 2 & 2 & 5 & 5 & 2 & 2 & 2 & 2 \\" ///
"Mean  Dep. Var. Control  & ${LossSit2Km_r_mean}  & ${LossSit2Km_r_mean}  & ${LossSit5Km_r_mean}  & ${LossSit5Km_r_mean}   & ${bosqprop2Km_mean}  & ${bosqprop2Km_mean} & ${bosqprop2KmC_mean} & ${bosqprop2KmC_mean} \\") 
eststo clear
************ Table B.5: Effects of eco-tourism promotion on other businesses' outcomes ******************************

preserve
set obs 10000
global var "M1_311 cost_wages M1_310 month_income_all"

* ANCOVA without controls 
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var'_nc: reghdfe `var' treatXtime `var'_m1, absorb(pair) cluster(codmpio)
qui sum `var'  if treat==0 & e(sample) & time == 1
glo  `var' _mean: di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA_nc: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo ancova_`var': reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM  DdistanciaBogm homicidios, absorb(pair) cluster(codmpio)
qui sum `var' if treat==0 & e(sample) & time == 1
glo  `var'_mean : di %6.2f r(mean)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_ANCOVA: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT without controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat if time == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var'_nc: reghdfe `var' treatXtime if time == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT_nc: di %6.2f `r(N)'/`Nrandos'
}
restore

* RCT with controls 
preserve
set obs 10000
foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment
reghdfe `var' randotreat edu_anos distanciaBogNM  DdistanciaBogm homicidios if time == 1, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind
}

*I do the actual regression
eststo rct_`var': reghdfe `var' treatXtime edu_anos distanciaBogNM  DdistanciaBogm homicidios if time == 1, absorb(pair) cluster(codmpio)

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))

*Randomization inference p-value two-tails
glo rand_`var'_RCT: di %6.2f `r(N)'/`Nrandos'
}
restore

* ANCOVA:
* output: Tables/others_outcomes 
estout ancova_M1_311_nc ancova_M1_311 ancova_cost_wages_nc ancova_cost_wages ancova_M1_310_nc ancova_M1_310 using "$outcomes/others_outcomes_business.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a, fmt(%12.2f ) labels ( "Adjusted R$^2$") )  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") prefoot("\bottomrule" ///
"RI p-value & ${rand_M1_311_ANCOVA_nc} & ${rand_M1_311_ANCOVA}  & ${rand_cost_wages_ANCOVA_nc}  & ${rand_cost_wages_ANCOVA}  & ${rand_M1_310_ANCOVA_nc}  & ${rand_M1_310_ANCOVA} \\" )

* RCT:  
* output: Tables/others_outcomes_rct
estout rct_M1_311_nc rct_M1_311 rct_cost_wages_nc rct_cost_wages rct_M1_310_nc rct_M1_310 using "$outcomes/others_outcomes_business_rct.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats( N r2_a , fmt(%11.2gc %12.2f  ) labels ("Observations" "Adjusted R$^2$"))  ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none) keep(treatXtime) varlabels( treatXtime "Promoted") ///
prefoot("\bottomrule" ///
"RI p-value & ${rand_M1_311_RCT_nc} & ${rand_M1_311_RCT}  & ${rand_cost_wages_RCT_nc}  & ${rand_cost_wages_RCT} & ${rand_M1_310_RCT_nc}  & ${rand_M1_310_RCT}  \\" ///
"Controls & No & Yes & No & Yes & No & Yes \\" ///
"Mean  Dep. Var. Control  & ${M1_311_mean }  & ${M1_311_mean }  & ${cost_wages_mean } & ${cost_wages_mean } & ${M1_310_mean } & ${M1_310_mean } \\")
eststo clear

****** Table B.6 Effects of eco-tourism promotion on businesses' and households' outcomes without controls  *******
* Regs from Table A.5

************* Table B.7: Effects of eco-tourism promotion on municipality deforestation
* output: Tables//defo_muni_reg.tex
preserve
bys codmpio time: gen num = _n
keep if num == 1

replace defo_muni = defo_muni/forest_2000*100
* ANCOVA 
eststo ancova_1: reghdfe defo_muni treatXtime defo_muni_m1, absorb(pair) cluster(codmpio)
qui sum defo_muni if treat==0  & time == 1
glo  mean_1: di %6.2f r(mean)

* ANCOVA with controls
eststo ancova_2: reghdfe defo_muni treatXtime defo_muni_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios, absorb(pair) cluster(codmpio)
qui sum defo_muni if treat==0 & time == 1
glo  mean_2: di %6.2f r(mean)

* RCT
eststo RCT_1: reghdfe defo_muni treatXtime if time == 1, absorb(pair) cluster(codmpio)
qui sum defo_muni if treat==0 & time == 1
glo  mean_3: di %6.2f r(mean)

* RCT with controls
eststo RCT_2: reghdfe defo_muni treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1, absorb(pair) cluster(codmpio)
qui sum defo_muni if treat==0& time == 1
glo  mean_4: di %6.2f r(mean)

* DiD
eststo DiD_1: reghdfe defo_muni treatXtime, absorb(codmpio time) vce( cluster (codmpio##time NIM)) 
qui sum defo_muni if treat==0 & time == 1
glo  mean_5: di %6.2f r(mean)
                                                                                               
* DiD with controls
eststo DiD_2: reghdfe defo_muni treatXtime edu_anosXtime distanciaBogNMXtime DdistanciaBogm homicidiosXtime, absorb(codmpio time) vce( cluster (codmpio##time NIM)) 
qui sum defo_muni if treat==0 & time == 1
glo  mean_6: di %6.2f r(mean)

restore              
estout  ancova_1 ancova_2 RCT_1 RCT_2 DiD_1 DiD_2 using "$outcomes/defo_muni_reg.tex", label replace style(tex)  starl(* 0.10 ** 0.05 *** 0.01)  stats(r2_a  N  , fmt(%12.2f  %11.2gc ) labels ("Adjusted R$^2$" "Observations" )) keep(treatXtime) ///
cells(b(star fmt(a2)) se(par fmt(a2))) mlabels(none) collabels(none)  varlabels(treatXtime "After X Promoted") ///
prefoot("\bottomrule" ///
"Controls & No & Yes & No & Yes & No & Yes \\" ///
"Mean  Dep. Var. Control  & ${mean_1}  & ${mean_2}  & ${mean_3}    & ${mean_4}  & ${mean_5}  & ${mean_6} \\") 
eststo clear

********************* FIGURES ************************
****** Figure A.1: Location of treatment and control municipalites ******
* Images/MapaTreatControl.png

****** Figure A.2: Timeline ******
* Images/Timeline.png

****** Figure A.3: Example of tourism products ******
* Images/PaginaAwake_lq.png

****** Figure A.4: Inference randomization ******
* Images/ntourist_histrando.pdf
preserve
set obs 10000
global var "ntourist"

foreach var of global var {
capture drop randonum  terand randotreat 
*Gen a variable where to store the treatment effects of the different randomization
gen terand=.
*Gen a variable where to store the treatment assignment of each randomization
gen randotreat=.
*Gen a variable to store the randon number that will guide the randomization
gen randonum=.
*decide how many randomizations
local Nrandos=10000
*set obs `Nrandos'
*Start the loop of randomizations
 qui forval i =1/`Nrandos' {
 *Set the seed on each iteration to get different randomization
set seed `i'
*I gen an indicator by group (the strata of randomization) to indicate if they 
*switch or maintain original treatment status
*The conditional guarantees I generate one per Group
sort pair
replace randonum=(runiform()<0.5) if pair!=pair[_n+1]
*I copy thi for the rest of the group
bys pair: egen randoind=max(randonum)
*Mathematical trick to maintain or switch treatment status
replace randotreat=1-(treat*randoind+(1-treat)*(1-randoind))

*Estimate the reg with the randomize treatment


reghdfe `var' randotreat `var'_m1  edu_anos distanciaBogNM DdistanciaBogm homicidios, absorb(pair)   cluster(codmpio)

*Store the ``treatment effect'' of the randomization
replace terand=_b[randotreat]/_se[randotreat] in `i'
drop randoind

}

*I do the actual regression
reghdfe `var' treat `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios , absorb(pair) cluster(codmpio)
*reghdfe `var' awake   $cond & time==1, absorb(Grupo)   cluster(codmpio)

*I count in how many cases I obtain a large estimate
count if terand>_b[treat]/_se[treat] & terand!=.

*Randomization inference p-value one tail
disp `r(N)'/`Nrandos'
texdoc init "$outcomes/ntourist_RI_pvalue_onesided.tex" , replace force 
local oneside: disp %4.2fc `r(N)'/`Nrandos'
tex `oneside'
texdoc close

*I count in how many cases I obtain a more adverse estimate
count if ((terand>abs(_b[treat]/_se[treat]) & terand!=.) | (terand<-abs(_b[treat]/_se[treat]) & terand!=.))
local line: di  %4.3f _b[treat]/_se[treat]

*Randomization inference p-value two tails
disp `r(N)'/`Nrandos'
texdoc init "$outcomes/ntourist_RI_pvalue_twosided.tex" , replace force 
local twoside: disp %4.2fc `r(N)'/`Nrandos'
tex `twoside'
texdoc close

twoway( hist terand, bin(30) fraction color(midgreen*0.4) fcolor(midgreen*0.5)), ///
       xline(`line', lc(red*0.7) lpattern(solid) lwidth(medium)) graphregion(color(white) ) xlabel(-6(2)6) plotregion(color(white) margin( b=0)) ///
xtitle("Randomization t-student", size(large))  ytitle("Fraction", size(large )) title("") legend(off)

graph export "${graphs}/`var'_histrando.pdf", replace
*graph export "${graphs2}/`var'_histrando.pdf", replace
}
restore

****** Figure A.5: Deforestation  effect by distance donus ******

* 1. Eco-tourism sites
preserve

local counter = 0
foreach var of varlist LossSitio1Km LossSitio1Km_2Km LossSitio2Km_3Km LossSitio3Km_4Km LossSitio4Km_5Km LossSitio5Km_6Km LossSitio6Km_7Km LossSitio7Km_8Km LossSitio8Km_9Km LossSitio9Km_10Km {
	
	* ANCOVA
	reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)
	estimates store reg_`counter'_1
	summ `var' if treatXtime  == 0 & e(sample)
	glo mean_`counter'_1 = `r(mean)'
	
	
	
	* RCT
	reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_LossSitio2Km == 1, absorb(pair) cluster(codmpio)
	estimates store reg_`counter'_2
	summ `var' if treatXtime  == 0 & e(sample) & time == 1
	glo mean_`counter'_2 = `r(mean)'
	
	local counter = `counter' + 1
	
}


coefplot reg_*, vertical gen(_c) name(coefplot1, replace) keep(treatXtime)

decode _cplot, gen(_ctemp)
gen _cnum_1tmp = substr(_ctemp, 5, 1)
gen _cnum_1 = real(_cnum_1tmp)
drop _cnum_1tmp
gen _cnum_2 = substr(_ctemp, 7, 1)

replace _cnum_1 = _cnum_1 + 0.2 if _cnum_2 == "2"

gen _coef_as_pct_mean_1 = 0
	foreach num in 0 1 2 3 4 5 6 7 8 9 {
	replace _coef_as_pct_mean_1 = _cb/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	replace _cul1 = _cul1/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	replace _cll1 = _cll1/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	}
	
gen _coef_as_pct_mean_2 = 0
	foreach num in 0 1 2 3 4 5 6 7 8 9 {
	replace _coef_as_pct_mean_2 = _cb/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	replace _cul1 = _cul1/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	replace _cll1 = _cll1/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	}
	
scatter _coef_as_pct_mean_1 _cnum_1 if _cnum_2 == "1", mcolor(navy) msymbol(S)  /// (Ready)
	|| scatter _coef_as_pct_mean_2 _cnum_1 if _cnum_2 == "2", mcolor(maroon)  ///
	|| rcap _cul1 _cll1 _cnum_1 if _cnum_2 == "1",  lc(navy) ///
    || rcap _cul1 _cll1 _cnum_1 if _cnum_2 == "2",  lc(maroon) /// 
    ||, xtitle(Bandwidth) xlabel(0 "0-1km" 1 "1-2km" 2 "2-3km" 3 "3-4km" 4 "4-5km" 5 "5-6km" 6 "6-7km" 7 "7-8km" 8 "8-9km" 9 "9-10km" ) ylabel(-200(50)200) ///
	graphregion(color(white)) graphregion(color(white)) xtitle("Donut distance (km)") ytitle("Effect size (%)") legend(order(1 "ANCOVA"  2 "Diff-in-mean"))
	
	graph export "$graphs/Plot_effectsize_donut_LossSitio.pdf", as(pdf) replace  
restore 

* 2. Business/HH
preserve

local counter = 0
foreach var of varlist LossCasa1Km LossCasa1Km_2Km LossCasa2Km_3Km LossCasa3Km_4Km LossCasa4Km_5Km LossCasa5Km_6Km LossCasa6Km_7Km LossCasa7Km_8Km LossCasa8Km_9Km LossCasa9Km_10Km {
	
	* ANCOVA
	reghdfe `var' treatXtime `var'_m1 edu_anos distanciaBogNM DdistanciaBogm homicidios if in_sam_LossCasa2Km == 1, absorb(pair) cluster(codmpio)
	estimates store reg_`counter'_1
	summ `var' if treatXtime  == 0 & e(sample)
	glo mean_`counter'_1 = `r(mean)'
	
	
	
	* RCT
	reghdfe `var' treatXtime edu_anos distanciaBogNM DdistanciaBogm homicidios if time == 1 & in_sam_LossCasa2Km == 1, absorb(pair) cluster(codmpio)
	estimates store reg_`counter'_2
	summ `var' if treatXtime  == 0 & e(sample) & time == 1
	glo mean_`counter'_2 = `r(mean)'
	
	local counter = `counter' + 1

}

coefplot reg_*, vertical gen(_c) name(coefplot1, replace) keep(treatXtime)

decode _cplot, gen(_ctemp)
gen _cnum_1tmp = substr(_ctemp, 5, 1)
gen _cnum_1 = real(_cnum_1tmp)
drop _cnum_1tmp
gen _cnum_2 = substr(_ctemp, 7, 1)

replace _cnum_1 = _cnum_1 + 0.2 if _cnum_2 == "2"

gen _coef_as_pct_mean_1 = 0
	foreach num in 0 1 2 3 4 5 6 7 8 9 {
	replace _coef_as_pct_mean_1 = _cb/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	replace _cul1 = _cul1/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	replace _cll1 = _cll1/${mean_`num'_1}*100 if _ctemp == "reg_`num'_1"
	}
	
gen _coef_as_pct_mean_2 = 0
	foreach num in 0 1 2 3 4 5 6 7 8 9 {
	replace _coef_as_pct_mean_2 = _cb/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	replace _cul1 = _cul1/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	replace _cll1 = _cll1/${mean_`num'_2}*100 if _ctemp == "reg_`num'_2"
	}
	
scatter _coef_as_pct_mean_1 _cnum_1 if _cnum_2 == "1", mcolor(navy) msymbol(S)   /// 
	|| scatter _coef_as_pct_mean_2 _cnum_1 if _cnum_2 == "2", mcolor(maroon) ///
	|| rcap _cul1 _cll1 _cnum_1 if _cnum_2 == "1",  lc(navy)  ///
    || rcap _cul1 _cll1 _cnum_1 if _cnum_2 == "2",  lc(maroon) /// 
    ||, xtitle(Bandwidth) xlabel(0 "0-1km" 1 "1-2km" 2 "2-3km" 3 "3-4km" 4 "4-5km" 5 "5-6km" 6 "6-7km" 7 "7-8km" 8 "8-9km" 9 "9-10km" ) ylabel(-200(50)200)  ///
	graphregion(color(white)) graphregion(color(white)) xtitle("Donut distance (km)") ytitle("Effect size (%)") legend(order(1 "ANCOVA"  2 "Diff-in-mean"))
	graph export "$graphs/Plot_effectsize_donut_LossCasa.pdf", replace  
restore 

****** Figure A.6: Histograms nworkers ******
twoway (histogram nworkers if in_sample_busoutcomes == 1 & treat == 0, percent width(2) start(0) color(red%30)) ///        
       (histogram nworkers if in_sample_busoutcomes == 1 & treat == 1, percent width(2) start(0) color(blue%30)), ///   
       legend(order(1 "Control" 2 "Treatment") rows(1)) xtitle("N workers") ytitle("%") graphregion(fcolor(white) lcolor(white))
	   graph export "${graphs}/Histogram_nworkers.png", replace
	   
***** Figure B.1: Distribution of crops and number of tourist ***** 
preserve 
use "${base_out}/Temporary/tourism_crop_by_month.dta", clear
twoway  (connected pct month,  lcolor(orange*0.5)   msymbol(none)  lwidth(medthick) lpattern(solid)) 				    ///
		, xlabel(1(1)12) xtitle("Month") ytitle("Percentage of crops") graphregion(fcolor(white) lcolor(white) margin(r=10 l=2 t=3 b=3))  				///
        ylabel(0(10)80, angle(h) nogrid)  	///
		graphregion(fcolor(white) margin(r=5 l=2)) plotregion(fcolor(white) lcolor(white))  legend(label(1 "Crops") rows(1))             
	    graph export "$graphs/graph_crops_pct.pdf", replace               

twoway  (connected tourism month,  lcolor(navy*0.5)   msymbol(none)  lwidth(medthick) lpattern(solid)) 				    ///
		, xlabel(1(1)12) xtitle("Month") ytitle("Percentage of companies") graphregion(fcolor(white) lcolor(white) margin(r=10 l=2 t=3 b=3))  				///
        ylabel(0(10)80, angle(h) nogrid)  	///
		graphregion(fcolor(white) margin(r=5 l=2)) plotregion(fcolor(white) lcolor(white))  legend(label(1 "Tourism") rows(1))             
	    graph export "$graphs/graph_tourism_pct.pdf", replace
restore

****** Figure B.2: Deforestation across years  ******
* Images/graph_defo.pdf 

preserve 
use "$base_out/defo_muni.dta" , clear
rename M_code codmpio
merge m:1 codmpio using  "$base_out/Temporary/treatment.dta", nogen keep(1 3)
replace treat=2 if treat==.
collapse (mean) def_rate, by(treat year)
label var def_rate "Deforestation rate"

twoway  (connected def_rate year if treat==0,  lcolor(green*0.5)   msymbol(none)  lwidth(medthick) lpattern(dash)) 				    ///
        (connected def_rate year if treat==1,  lcolor(green) msymbol(none) lwidth(medthick) lpattern(solid))					    ///
        (connected def_rate year if treat==2,  lcolor(green*1.5)  msymbol(none) lwidth(medthick) lpattern(dot)) 					///
		, xlabel(2000(2)2018) xtitle("")  graphregion(margin(r=10 l=2 t=3 b=3))  				///
	    legend(on stack region(lcolor(white)) size(small) cols(4) label(1 Control) label(2 Treatment) label(3 Rest of country) rows(1)    ///
		  symx(23) symy(2))  ylabel(0(0.1)0.6, labsize(vsmall) angle(h) nogrid)  	///
		graphregion(fcolor(white) margin(r=5 l=2)) plotregion(fcolor(white) margin(r=6 l=0 t=0 b=0) lcolor(white))                  
	    graph export "$graphs/graph_defo.pdf", replace               
restore

****** Figure B.3: Buffer ******
*Images/buffer.pdf

****** Figure B.4: Buffers with forest raster ******
*Images/bufferb.PNG

****** Figure B.5: Buffers for calculating externalities ******
*Images/spillovers_map.pdf

****** Figure B.6: Deforestation balance years before treatment ******

preserve 
use "${base_out}/balance_data_for_reg.dta", clear

* to check outliers
* summ LossCasa2Km if year == `n' & in_sam_LossCasa2Km == 1, d
* summ LossSitio2Km if year == `n' & in_sam_LossCasaSitio

local counter = 0
foreach n of numlist 2001(1)2017 {
	reghdfe LossCasa2Km treat distanciaBogNM  DdistanciaBogm homicidios if year == `n' & in_sam_LossCasa2Km == 1, abs(pair) cluster(codmpio)
	estimates store reg_`counter'_1
	reghdfe LossSitio2Km treat distanciaBogNM  DdistanciaBogm homicidios if year == `n' & in_sam_LossSitio2Km == 1, abs(pair) cluster(codmpio)
	estimates store reg_`counter'_2
	
	local counter = `counter' + 1
}

coefplot reg_*, vertical gen(_c) name(coefplot1, replace) keep(treat)

gen num_1 = mod(_cplot, 2)
gen num_2 = ceil(_cplot/2)-1
replace num_2 = num_2 + 0.2 if num_1 == 1

scatter _cb num_2 if num_1 == 1, mcolor(navy)  /// 
	|| scatter _cb num_2 if num_1 == 0, mcolor(maroon) ///
	|| rcap _cul1 _cll1 num_2 if num_1 == 1,  lc(navy)  ///
    || rcap _cul1 _cll1 num_2 if num_1 == 0,  lc(maroon) /// 
    ||, xtitle(Bandwidth) xlabel(0 "2001" 1 " " 2 "2003" 3 " " 4 "2005" 5 " " 6 "2007" 7 " " 8 "2009" 9 " " ///
	10 "2011" 11 " " 12 "2013" 13 " " 14 "2015" 15 " " 16 "2017") ylabel(-1(0.25)1) ///
	graphregion(color(white)) graphregion(color(white)) xtitle("Deforestation (Ha)") ytitle("Deforestation (Ha)") legend(order(1 "Business/HH"  2 "Eco-tourism sites"))
	graph export "$graphs/Balance_Loss.pdf", replace  
estimates clear 
restore 




  